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In this paper we analyze the predictions of the forward approximation in some models which 
exhibit an Anderson (single-) or many-body localized phase. This approximation, which consists 
in summing over the amplitudes of only the shortest paths in the locator expansion, is known to 
over-estimate the critical value of the disorder which determines the onset of the localized phase. 
Nevertheless, the results provided by the approximation become more and more accurate as the local 
coordination (dimensionality) of the graph, defined by the hopping matrix, is made larger. In this 
sense, the forward approximation can be regarded as a mean field theory for the Anderson transition 
in infinite dimensions. The sum can be efficiently computed using transfer matrix techniques, and 
the results are compared with the most precise exact diagonalization results available. 

For the Anderson problem, we find a critical value of the disorder which is 0.9% off the most 
precise available numerical value already in 5 spatial dimensions, while for the many-body localized 
phase of the Heisenberg model with random fields the critical disorder he = 4.0 ± 0.3 is strikingly 
close to the most recent results obtained by exact diagonalization. In both cases we obtain a critical 
exponent i/ — 1. In the Anderson case, the latter does not show dependence on the dimensionality, 
as it is common within mean field approximations. 

We discuss the relevance of the correlations between the shortest paths for both the single- and 
many-body problems, and comment on the connections of our results with the problem of directed 
polymers in random medium. 


I. INTRODUCTION 

The propagation of waves and quantum particles in a 
disordered medium is a fascinating and challenging prob¬ 
lem in statistical mechanics, with plenty of relevance 
for experimentJi^. Among the phenomena that occur 
in such a set ting , Anderson localization is probably the 
most strikinglEl. “Anderson’s theorem” in Ref. [3] states 
that for sufficiently strong disorder, diffusive transport 
is completely suppressed in single particle problems on 
a lattice. The study of the transition that separates the 
two phases (localized and delocalized) has resisted an ex¬ 
act solution for about 60 years, and numerical methods 
are still at the core of the advances in this topieP. 

Recently, interest in Anderson localization has surged 
due to the work by Basko, Aleiner and Altshuleil^, hence¬ 
forth denoted by BAA. There, the phenomenon of many- 
body localization (MBL), i.e. the stability of the Ander¬ 
son insulator to the addition of interactions, is investi¬ 
gated. BAA’s work has been extended and re-interpreted 
in several other work^^HIIl^ ^nd MBL appears now to be 
the most robust mechanism to break the ergodicity that 
is typical of generic interacting systems. 

The core of BAA’s analysis relies on perturbatively ac¬ 
counting for the interactions, at finite temperature and 
particle density. Technically, they consider the pertur¬ 
bation theory for the imaginary part of the propagator 
of an excitation on top of an eigenstate by means of the 


Keldysh formalism. For sufficiently weak interactions, 
the perturbative series is shown to converge with prob¬ 
ability equal to one. In the spirit of Ref. H] this implies 
the localization of the excitations themselves, and the 
absence of transport. 

As an alternative route, the MBL problem can be in¬ 
terpreted as a single particle tight-binding problem in 
the space of many-body configurationJi^, with the inter¬ 
actions playing the role of an effective hopping. However, 
several issues arise in this formulation. First, the on-site 
energies in the resulting effective lattice are not inde¬ 
pendent variables drawn from the same distribution, but 
they are strongly correlated. Secondly, the connectivity 
of a configuration in the many-body problem scales with 
(a power of) the system size: it diverges in the thermody¬ 
namic limit, and thus it is impossible to define a limiting 
graph. As a consequence, one needs to define an effective 
connectivity which stays of 0(1) in the thermodynamic 
limit (a similar phenomenon is observed in Ref. mi. Fi¬ 
nally, the number of paths of a given length connecting 
two many-body configurations grows factorially with the 
distance between them, the distance being defined as the 
minimum number of actions of the interaction operator 
needed to connect the initial to the final configuration. 
Since the distance between two configurations can be of 
the order of the system size, the number of paths can 
grow factorially in the system size. When the mapping 
to a single particle problem is applied to a system of 


2 


N interacting spins, it results into a correlated-disorder 
problem on a (section of) an iV-dimensional hypercube. 

De spite these complications, in some recent 
workJSMIill this approach has been used success¬ 
fully to estimate, among other things, the boundaries of 
the MBL region. The analytical calculations in Refs. [TH 
and HU were shown to be in very good agreement with 
the numerical results in the same papers, obtained 
with exact diagonalization. The analytic results are 
derived revisiting an approximation already used in 
Refs. 11811201 and HU which consists in calculating the 
Green’s functions by retaining only the lowest order in 
the hopping. Recently, this approximation has been used 
in Ref. [22] to derive the power-law tail of the distribution 
of the wave function amplitudes on a Bethe lattice in 
the localized phase. 

In this paper we discuss in detail this approximation, 
dubbed “forward approximation” (FA), by illustrating 
its virtues and limitations and its connections to other, 
seemingly unrelated, problems of statistical physics. The 
paper is organized as follows: in Section [n] we derive the 
expression for the wave function amplitudes in forward 
approximation, and discuss a criterion for localization 
given in terms of the probability of resonances. In Sec- 
tion jlll] after recalling some known results on the Ander¬ 
son model on the Bethe lattice, we focus on d-dimensional 
systems (with d = 3—6). We compare the analytic results 
we get for the critical disorder within the FA with the 
numerics, showing how the approximation gives better 
results as the dimensionality d is increased. We then dis¬ 
cuss the application of the aforementioned technique to a 
many-body problem of interacting spins in a disordered 
environment. In Section IIVI we discuss the relevance of 
the correlations and interference between different paths 
connecting two configurations, both for the single parti¬ 
cle and for the many body problem. In the Conclusions 
we comment on the various possible directions in which 
this work can be extended, focusing in particular on the 
application of the forward approximation to the study 
of the MBL phase and of the many-body localization- 
delocalization transition. 


II. DERIVATION 

A. The forward approximation for the 
eigenfunctions 

To begin with, we derive the expression for the wave 
function amplitudes in FA for a single particle hopping 
on a finite lattice with on-site disorder. We consider the 
Hamiltonian: 


N 

iL = ^eic^c,(c^c^ +h.c.) . (I) 

(ij> 

In the Anderson model, the are independent random 
variables uniformly distributed in [—IF/2, IT/2], one for 


each of the N sites in the lattice. The edges {i,j) define 
the lattice geometry. The lattice constant is set to one, 
and we denote with L the length scale characterizing the 
size of the lattice (for a cubic lattice in dimension d the di¬ 
ameter is L = for a Bethe lattice or regular random 

graph it is L = In [{N — 1){K — 1)/{K -|- 1) -|- 1] / In AT « 
In N/ In K, where AT-I-1 is the connectivity of the lattice). 
The distance d{a, b) between two arbitrary sites a, b of the 
lattice is the number of edges in a shortest path connect¬ 
ing them. We refer to it as the lattice distance in the 
following. 

We consider the matrix elements of the resolvent be¬ 
tween two states associated to the sites a, b in the lattice: 

G{b,a,E) = {b\^^\a), (2) 

at energy E. They have the following expansion in series 
of t 

C,i,a.E) = T^^ E (3) 

pGpaths(a,6) i£p 


where the sum runs over all the paths p in the lattice 
connecting the sites a and b, and i labels the sites visited 
by the path p (site a excluded). Given two sites a,b at 
lattice distance n, the lowest orders in the expansion are 


G{b, a, E) 


lit t t 

E - ea E - ei E - 62 "' E - e„_i E - eb~^ 
+ other paths of length n + ■ ■ ■ 

( 4 ) 


where the sites (a, I, 2,..., n — 1,6) belong to one of the 
shortest paths connecting a and b. 

On the other hand, from the spectral decomposition it 
follows 


G{b, a, E) 


E 


E - Ea 


( 5 ) 


and thus the residue at E = E^ gives 


lim {E-Ea)G{b,a,E) =ilja{b)'ilJa{a), (6) 


assuming no degeneracy of the eigenvalues. 

The path representation ^ does not have a pole at E ^, 
to no order in t. To get the exact poles it is necessary 
to re-sum the closed paths in the series expansion. Once 
this is done, the full series is re-cast into a sum over the 
non-repeating paths paths* (a, 6) connecting the sites a 
and b: 


G{b, a, E) 


1 

E — ea ~ Fla(A') 

E 

pGpaths* (a,6) 



t 

-4^\e)' 


( 7 ) 


Here T,a{E) is the local self-energy at the site a, defined 
through the identity: 


G(a, a, E) = 


1 

E — ta — T,a{E) 


( 8 ) 
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It is equal to the sum of the amplitudes of all the closed 
paths in which site a appears only as starting and ending 
point; to lowest order in t 

^a{E)=Y, (9) 

T a Fin 


where da is the set of nearest neighboring sites of a. The 
path-dependent term is a modified self-energy, 

which re-sums the loops around site i, never crossing site 
i again, nor any of the sites (a, 1, • • • , z—1) already visited 
by the non-repeating path p. 

The expansion Q in non-repeating paths has several 
advantages. First, while paths{a, b) is an infinite set, 
even for a finite lattice, paths* {a, b) is finite for a finite 
lattice. Moreover, 0 is free of the divergences affecting 
([^ that are due to local resonances, i.e. to sites i,j at 
bounded distance satisfying |ei — ej| ~ A, w ith A small 
(we clarify the exact meaning of small in Sec. im Local 
resonances necessarily occur also in the localized phase, 
and produce arbitrarily large factors in ([^, correspond¬ 
ing to the paths repeatedly hitting the resonant sites an 
arbitrary number of times. These repetitions of large 
factors lead to the divergence of the naive perturbation 
series in t, but are re-summed into self energy corrections 
in Q. The “renormalized” expansion in non-rraeating 
paths is found to converge in the localized phas™, when 
resonances do not proliferate at asymptotically large dis¬ 
tances in space, and the hopping hybridizes only the de¬ 
grees of freedom in a finite, albeit possibly big, region of 
space. An analogous resummation procedure is discussed 
in Ref. 1161 where the perturbation theory for quasi-local 
conserved operators is shown to converge in the MBL 
phase. 

The expression for the eigenfunction is obtained from 
the resolvent as follows: the eigenenergy satisfies 
Ea = Ca + E.a{Ea), thus the first factor of Q has a 
pole at Ea with residue |z/^Q(a)P, as it follows from ^ 
and Q. Then: 


lim {E — Ea)G{b,a, E) = 

E—>-Ecc 

|'0a(a)P lim W TT - - — 

pGpaths*(a,f)) JGp I 


which gives 

i’aib) = Ipaia) 


E n 


pGpaths*(a,b) iGp Cz 


( 11 ) 

with z/>a(a) obtained from Ea{Ea) using Q.122I 

From © we can read the expression of the wave 
function amplitudes to lowest order in t. Assume that 
a labels an eigenstate localized at site a for t —>• 0. 
Since Sa = 0{F), we have to lowest order E^ —t Ca, 
ipaia) 0, giving 


M)= E (12) 

pGspaths(a,6) iGp 


where the set spaths(a, 6) C paths* (a, 6) contains the 
shortest paths from a to b. Note that this derivation 
does not rely on the particular structure of the lattice 
nor on the independence of the on-site energies, and thus 
it can be straightforwardly generalized to hopping prob¬ 
lems on graphs with different geometry or more general 
distribution of the random energies, such as the ones in 
Sec. lmCl 

The effect of the modified self energy corrections 
is to weaken the role of resonances. Indeed, let 
us assume that in the path p there is a site i which is 
resonant with a, jca — e^l A. In this case the forward 
approximation (12) will contain the very large term t/A. 
However, the correction I]-((.\(ea) in the previous site also 
contains such large term, leading to a compensation. As 
we shall see in the following, neglecting this effect leads to 
an overestimate of the minimum disorder strength needed 
to localize the system. On the other hand, the loops con¬ 
tributing to the self energy corrections become less rel¬ 
evant when the dimensionality (or connectivity) of the 
lattice is increased. Thus, the FA is expe cted t o give 
faithful results in higher dimension (see Sec. IIIB). 


B. Probability of resonances and criterion for 
localization 


For a single particle problem on a finite dimensional 
lattice, we define an eigenstate ipa of a system of size L 
localized if (with probability I over the disorder realiza¬ 
tions) the probability of finding a particle at a distance 
0{L) from the localization center of the state tends to 
zero in the limit of large L. More precisely, let a denote 
the localization center of ip a, and 


ipr = max \ipa(b)\. 

b: d{b,a)—r 


(13) 


We define the state ipa localized if there exists a finite 
^ > 0 such that 


V 2r - 2C 


—I for 


(14) 


Namely, we require that the random numbers \ipa(b)\'^ 
can be enclosed in an exponential envelope for all b suf¬ 
ficiently far from the localization center of the state. By 
means of the Kubo formula, it is possible to show that 
this condition on the eigenstates implies the vanishing of 
the DC conductivity. We identify the localization length 
of Ipa with the minimum value of ^ for which Eq. (14) 


is true. It does, in general, depend on the state; how¬ 
ever, it is supposed to depend smoothly on the energy 
Ea in the thermodynamic limit. A mobility edge exists 
whenever there is band of energies for which such mini¬ 
mum is not finite. In particular, at fixed energy and at 
the corresponding critical value of disorder W = Wc, the 
localization length diverges and the asymptotic bound 
in Eq. (14) ceases to hold for any finite This entails 
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that for any arbitrarily small, positive e = and at 
arbitrarily large distances r from the localization center, 
there exist sites b such that the ratio log |'0ct(&)P/2r’ ex¬ 
ceeds the constant — e with some finite probability. We 
expect the delocalized phase W < Wc to be characterized 
by the stronger condition: 


the problem amenable to analytic calculations! ^^ l ^° ^ We 
briefly recall some results in the following. 

Let a be the root of the tree, and K the branching 
number (the connectivity is -1-1). Within the FA we get 
that the wave function at one particular point at distance 
L from the root is given by 


log li’r 

2 r 


> —e I —>• 1 for r —)• 00 , 


(15) 


holding for any arbitrarily small, strictly positive value 
of e. 


From the equations (141, (15) it follows that the 


localization-delocalization transition can be detected an¬ 
alyzing the statistics of the wave function amplitudes as 
a function of distance. In the following, we compute the 
wave function amplitudes in FA and determine numer¬ 
ically the probability in Eq. (HH), choosing e of the or¬ 
der of the numerical precision. We refer to the resulting 
probability as the “probability of resonances”. The ter¬ 
minology is motivated by the fact that within the FA, 
an amplitude of 0(1) at a site b at distance r from the 
localization center a corresponds to a resonance between 
the two sites. Indeed, the corresponding two sites prob¬ 
lem can be considered as a two-level system with reduced 
Hamiltonian 


h = 


0 hr 

hr A 


where A = ea — Cb, and 

hr = t ^ 


pGpaths* (a,&) iGp 


Ca-ei- 


(16) 


(17) 


where the products are taken over all sites in the path, 
excluding a, b. The sites are resonant when the energy 
difference Eq — et is small, and precisely |A| < hr- Con¬ 
sidering hr to lowest order in t, one finds that this is 
equivalent to |^/’(5)| > 1, with ipib) computed in the low¬ 
est order FA. Thus, with (141 and (15) one is probing 


the statistics of resonances within the FA, and requiring 
that the probability to find at least a resonant site at any 
sufficiently large distance r from the localization center 
decays to zero in the localized phase. This criterion in¬ 
volving resonances allows to obtain an estimate for the 
critical disorder within the FA also on the Bethe lattice. 


where the exact eigenstates satisfy Eq. (14) even in the 
delocalized phase, due to normalization. 


III. RESULTS 

A. Warming up: the Bethe lattice 

The simplest setting for the application of the forward 
approximation is a Bethe lattice: in this case, given two 
sites a, b there is only one non-repeating path connecting 
them, along which all the energies are i.i.d. . This makes 


^L = Yly^. (18) 

i—l 

This random variable has a power-law tail distribution 
for any distribution of Ci having a support S such that 
Ea G S, as one can see from the divergence of the first 
moment of the absolute value of ^pL■ It is convenient to 
consider the distribution of the logarithm oi ip l, whose 
moments are all finite. Let us remind the reader that 
we choose G [—W/2,W/2], and in this calculation for 
simplicity we set Cq = 0. Defining 

a;L = In =-2^1n(|ej|/t) (19) 


we find that 


{xl) = 2L\n{2et/W). 


( 20 ) 


The ratio {xl) /L is the typical decay of wave function 
amplitudes from the origin of localization 


= 2ln{2etlW). 


( 21 ) 


However, this is not the localization length ^ as defined in 
(14). The latter is indeed a uniform bound over the full 
set of ^ points at distance L from the localization 
center: it is determined by the decay rate of the maximal 
amplitude over sites in each shell at distance L. This is a 
point-to-set correlation function decay, which is familiar 
in the study of disordered systems on the Bethe lattice 
(on a regular lattice the point-to-set is substituted by 
the point-to-point, but with exponentially many shortest 
paths leading to the final point). 

The typical value of the maximal amplitude x^ among 
samplings is the largest solution of 


''L ■ 


K^Pixl) ~ 1 . 


( 22 ) 


Here P{x) is the distribution of xl, and the paths 
are treated as independent. For large L, what matters is 
the tail of the distribution P{x). If we rescale 


z=^+HW/2t), 


(23) 


we get, for large z: 


P{z) ~ exp {—L{z — 1 — In z)). 


(24) 


The probability distribution of z can be found inverting 
its Laplace transform, which is computable since z is a 
sum of i.i.d. variables. The maximum z* over sam¬ 
plings of z is the solution of 

0 = -z* + l + \nz* + \nK = ln(z*eAe-^*), (25) 
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FIG. 1. (Color online) Anderson model on a cube of side L = 
3. The red, dashed edges form one of the non-repeating paths 
connecting the sites a and b, of length r = {L — l)d = 6. The 
other elements in the set paths* (a, b) are obtained following 
the arrows. 


and ~ —2Lln(VF/2I) -|- 2Lz*. Localization occurs as 
long as x^/L < 0, so the critical condition can be written 
as 

ln(A:ez*e"^'') =0, (26) 

z* - ln(We/2t) = 0. (27) 

By eliminating z* we recover the familiail^ equation: 

Wc = 2teKln{Wj2t). (28) 


treat as the origin. Given any other site b, the orientation 
of the non-repeating, shortest paths from o to & induces 
a natural orientation of the edges of the cube, which is 
thus directed, see Fig.[^ Let r be the lattice distance of 
the sites with respect to the origin a: for a cube of side 
L, the maximum r is rmax = {L — l)d, corresponding to 
the site at the opposite corner of the cube with respect 
to a. 

In an infinite cube the number of points at distance r 
from the origin is (r-l-d—I)!/(((i—I)!r!), and thus it grows 
at most polynomially in r, slower than ~ r‘^. Naively, one 
would be led to think that the transition is given by the 
divergence of ^typ, see Eq. (21). However, the number 
of minimum length paths leading from the origin to an 
arbitrary site at lattice distance r scales exponentially 
with r, as ~ (F. Therefore, unlike in the Bethe lattice 
case, in finite dimension the wave function amplitude at 
a given site 6 is a sum over exponentially many correlated 
terms 


pGspaths(a,D) 

where 

pespaths(a,b) i—1 ^ ^ 


and the random variables are uniformly distributed in 
[—1/2,1/2]. In the following, we consider the probability 
distribution of the random variable 


Moreover, z* is given in terms of W/, as 

z* =ln(W',/2t), (29) 

so we can write the localization length as 

= -xl/L = 2MW/W,). (30) 

This gives a mean-field exponent for the divergence of 
the localization length at the transition: 


, ^ \0gWr? 
2r 


(34) 


for different values of r. Here denotes the maximum 
among all the rescaled amplitudes (33) at sites that are 
at lattice distance r with respect to the origin a. The 
probability of resonances for arbitrary values of t and 
W, see Eq. (15), is easily recovered from the cumulative 
distribution function of Zr as: 


2Wc 

^ “ jVE-lEer 



Note that ^typ < C irrespective of the value of the disor¬ 
der. The belief expressed in Ref. HHis that this behavior 
persists even within the delocalized regime, in the form 
of very irregular (multi fractal) eigenstates. Note addi¬ 
tionally that the difference between ^typ and ^ is due to 
the exponential sampling of the probability distribution, 
and it extends both to the finite-dimensional cube and 
to the many-body case. 


B. d-dimensional cube 

Consider now the case of a d-dimensional lattice of side 
L, and let a be the site at one corner of the cube, which we 


with e arbitrarily close to zero. According to Eqs. ([T4|) 
and 


( |15[ ), the density of Zr becomes asymptotically 
peaked at log(Wc/t) for r —oo, with width going to 
zero with r. Thus, the critical value of disorder can be 
estimated inspecting the scaling with r of the probability 
density of Zr. 

The distribution of Zr is hard to determine analyti¬ 
cally, due to the correlation between the different short¬ 
est paths. To account for such correlations, we compute 
the amplitude (33) numerically by means of a transfer 
matrix technique, and use the resulting values to deter¬ 
mine the probability (35) with e smaller than the numer¬ 
ical precision. The convenience of the transfer matrix 
method relies on the fact that it takes only polynomial 
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time in r, as it was realized by Medina and KardaJ^ in 
their treatment of the Nguyen, Spivak, and ShklovskiP^ 
(NSS) modelP 

The numerical computation is as follows. We fix t = 1 
and introduce the matrix T defined as 


T = WAf, (36) 

where is the forward adjacency matrix of the lattice 
(that is, the adjacency matrix associated to the directed 
cube of FigQ, and W is a diagonal matrix with compo¬ 
nents: 


W = diag 


- ef 


(37) 


k / k=l,...L'^ 


We initialize the system in the state = |a) com¬ 

pletely localized in the origin a, and iteratively apply the 
transfer matrix T. A single iteration gives 

^ l^i) + \h) +..., (38) 


where Zi, • • • ,ld are the forward neighbors of site a. The 
value of ipa{b) equals ipaib) = {bl’tp'-^^), where |6) is the 
state completely localized in the site b and r is the lattice 
distance between a and b. 


We fix = 0 and compute the rescaled amplitude (33) 


for all the points 6 on a shell at the same lattice distance 
r = Tmax — c from the origin of a hypercube of side L. 
Here c ~ 0(1) is fixed so as to have about 20 points per 
each size of the hypercube. We determine the maximal 
'0(. among the wave function amplitudes on those sites. 
We repeat the procedure for hypercubes of different sizes, 
with 0(10®) disorder realizations for most system sizes, 
decreasing to 0(10®) realizations only for the biggest sys¬ 
tem sizes that we consider (e.g. in d = 3 we take system 
sizes r = 10 through 292, with 1.5 • 10® disorder real¬ 
izations up to r = 202 and 2.5 • 10® realizations up to 
r = 292). 

As we discuss in Section [Tvl the main contribution to 
the transfer matrix result comes from only one of the ex¬ 
ponentially many paths in (33), and the results obtained 


with the transfer matrix technique are faithfully repro¬ 
duced by analyzing the statistics of the dominant path 
alone. The latter can be determined (see Section IV) with 
an algorithm that is computationally more efficient than 
the transfer matrix, allowing to access to much bigger 
system sizes. The results presented in this sections for 
d = 6, 7, as well as for the higher values of r in d = 3 — 5, 
are obtained with this procedure. 


1. Fluctuations of the wave function amplitudes 

In Fig. we plot the probability density of the variable 
Zr defined in Eq. ( |M| , for different values of r in d = 3. 
The plot shows a drift of the position of the peaks with 
increasing r, together with the shrinking of the width of 



FIG. 2. (Color online) Probability density of the variable Zr 
defined in Eq. ( |34[ ), for different r and d = 3. For r —> 00 , the 
curves become peaked around the critical value log(IFc/t). In¬ 
set: cumulative distribution function. Each curve is obtained 
with 1.5 • 10® disorder realizations. Very similar results are 
obtained for higher dimensionality. 


the distribution, in agreement with the conditions (141, 
(15). Plots of the r-dependence of the variance of 
(34) are given in Fig. in log-log scale for d = 3 — 6. 
The linear behavior indicates that the fluctuations of Zr 
decay to zero as a power law in r, with a coefficient that 
depends on the dimensionality. The higher cumulants of 
the distribution exhibit a similar linear behavior in log- 
log scale. Moreover, the numerical computation indicates 
that for fixed d the probability densities of the variable 


Zr - (Zr) , , 

Z,. = ^ (39) 

<JZr 

collapse to a limiting curve for increasing r, see Fig. 

As shown in the same plot, for fixed r and varying di¬ 
mensionality, the distribution of Zr does not change sig¬ 
nificantly, except for a weak d-dependence of the tails. 

These numerical observations are compatible with the 
following large r scaling form for Zr- 


rZr 


~ r log 



-b 


(40) 


where m is a random variable of 0(1) with a distribution 
which depends weakly on the dimensionality. 

According to (40), for large r the fluctuations a\ de¬ 
cay to zero with the power From the linear 

fit of log (o'z„) we extract the numerical estimate of the 
exponent in (40), which we denote with ujpA(d). The 
results are reported in Table 

In order to characterize the limiting distribution in 
Fig. 0 we compute the skewness Sk = K^j and the 
kurtosis Kur = ka! of the density of Zr (here Ki de¬ 
notes the i-th cumulant of the distribution). From ( |40| ) 
it follows that these parameters approach the ones cor¬ 
responding to the variable u in the limit of large r. We 
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FIG. 3. (Color online) Variance a\^ of the variable Zr de¬ 
fined in Eq. (34l. The plot is in log-log scale. The points 


corresponding to lar ger r are htted linearly, according to the 
scaling form Eq. (401, and the values of the exponents uiFA{d) 
reported in Table |I| are extracted from the coefficient of the 
linear term in the fit. The number of realizations is 1.5 • 10® 
for r smaller than 202, 53, 52, 40 for d = 3, 4, 5, and 6, re¬ 
spectively, and 2 ■ 10® for larger values of r. Inset Mean value 
of the variable rZr defined in Eq. (341, for d = 3. The fit 


is linear with a correction oc in agreement with the 

scaling form in Eq. ( |40[ ), with the value of ujfa{S) give n in 
Table The results of the fit are, with reference to Eq. (42l; 
Cl = -18.2 ± 0.3, Wc = 27.03 ± 0.02, C 2 = 29.6 ± 0.8. The 
same behavior holds for higher dimensionality and results in 
the estimates of the critical disorder values in Table m 


TABLE I. Values of the exponent ujFAid) governing the decay 
of the fluctuations of Zr with r, see Eq. ( |40[ ). A comparison 
is made with the values of the droplet exponents u)dp{D) 
obtained numerically for the directed polymer in dimension 
1 -I- (d — 1). The numerical values are taken from Appendix 
A in Ref. 1281 


d=D-tl 

UJFA{d) 

ludp{D) 

3 

0.278 ± 0.005 

0.244 

4 

0.23 ±0.01 

0.186 

5 

0.191 ±0.007 

0.153 

6 

0.168 ±0.006 

0.130 


restrict to d = 3, for which we have the largest statis¬ 
tics available. Plots of the r-dependence of Sk and Kur 
are given in Fig. The asymptotic values are estimated 
to be Sk = 0.34 ± 0.02 and Kur = 3.24 ± 0.04, see the 
caption of Fig. [5] for details. 


2. Estimate of the critical disorder 

To determine the critical value of disorder for t = 1, 
we extrapolate the asymptotic limit of the typical value 
of Zr- Since the distribution is not fat tailed, we can 
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FIG. 4. (Color online) Probability density P{Z) of the vari¬ 
able Zr dehned in Eq. (391. Each curve is obtained with 
1.5 • 10® realizations. Top. Density of Zr for different values 
r and d = 3. The curves seem to converge to a unique limit¬ 
ing distribution with increasing r. Bottom. Density of Zr for 
fixed r = 52 and different dimensionalities. 


equivalently consider the averages of Zr and set: 


(Zoo) = lim {Zr) = log {Wc). 

r—^oo 


(41) 


The inset in Fig.|^shows the scaling with r of r{Zr). The 
average grows linearly in r, in agreement with Eq. (401. 
We fit the data with the form 


{rZr) = Cl + \og{Wc) r + C2 , (42) 

with the numerical values uj{d) = ujpA{d) reported in 
Table |T). 

The resulting estimates of the critical disorder, which 
we denote with are displayed in Table [ll| For the 

smallest dimensions, a comparison is made with the criti¬ 
cal values VF™™ determined in Refs. and [501 bv means 
of a combination of exact diagonalization and transfer 
matrix techniques. 

The data in Table |ll] clearly show that the FA gives 
an upper bound to the critical disorder, since the renor¬ 
malization of the energy denominators provided by the 
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FIG. 5. (Color online) Skewness Sk of the distribution of the 
variable Zr defined in Eq. ( |39[ ), for d = 3. The red dashed 
line is a fit of the form a + /3r^, with a, /3, 7 free parameters. 
The coefficient a is the estimate of the asymptotic value of 
the skewness, and it equals Sk = 0.34T0.02. Inset. Kurtosis 
Kur of the distribution of Zr for d = 3, as a function of r. The 
fitting procedure is analogous to the one for the skewness, and 
results in Kur = 3.24 ± 0.04. 
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TABLE II. Comparison between the critical value for local¬ 
ization in the Anderson model in d dimensions predicted by 
the forward approximation and the numerical results 

of Ref. 1291 The relative error decreases faster than 
d“®, presumably exponentially. For d = 6 the transition value 
= 77.0±0.3 can be compared with the result of Ref. 1311 
= 74.5 ± 0.7. This number is however an underestima¬ 
tion of the transition due to the choice of boundary conditions. 
For 7 dimensions there is no available numerics to compare 
with. 


d 

IFf^ 

^num 

Error 

3 

27.03 ±0.03 

16.536 ±0.007 

39% 

4 

41.4 ± 0.1 

34.62 ± 0.03 

16% 

5 

57.8 ±0.2 

57.30 ±0.05 

0.9% 

6 

77.0 ±0.3 

- 

- 

7 

93.8 ±0.3 

- 

- 


(modified) self-energy corrections are neglected, and the 
effects of resonances are thus enhanced. However, in¬ 
creasing the dimensionality the discrepancy between the 
numerical estimates of Wc decreases; the enhanced pre¬ 
cision of the FA result is due to the fact that the loops 
giving rise to the self-energy corrections become less rel¬ 
evant in higher dimensional lattices, and thus the FA 
becomes asymptotically exact in this limit. 


3. Divergent length scales and critical exponents 


For fixed values of W and for finite r, the probabil¬ 


ity of resonances (351 is determined by the tails of the 
distribution of Zr- 

For increasing r, the asymptotic limit is approached 
in a different way at the two sides of the transition: for 
W > Wc, the probability of resonances goes to zero ex¬ 
ponentially with r. Below the transition, the probability 
converges to one much faster, with corrections that are 
only double exponential in r. We justify analytically this 
behavior in Appendix by computing an approximate 
expression for the density of the variable Zr- The ap¬ 
proximation consists in considering the different paths 
contributing to it as independent variables. 

Examples of the fits of the probability of resonances 
are shown in Fig. To extract a VF-dependent length 
scale 1{W), we perform an exponential fit of the form: 




1{W)\ 


(43) 


for W > Wc- For W < Wc we determine 1{W) by means 
of the linear fit 


log 


log 


logl^rP 

2r 



02 (IF) 


mr 

(44) 


The length scale 1{W) is plotted in Fig. for d = 3. 
We expect it to diverge in the same way as the lo¬ 
calization length/correlation length does in the local¬ 
ized/delocalized phase, respectively. We find that ^(IF) 
diverges as a power-law at a critical disorder compatible 
with the values of Wc^ listed in Table A fit of the 
form log {1{W)) = logc—izlog |IF—IF/’'^| results in an ex¬ 
ponent that is compatible with v ^ 1 for all dimensions, 
consistently with the Bethe lattice picture and with the 
results in Appendix]^ However, some deviations can be 
observed: a more careful analysis of the numerical data 
will be presented in a future publication. 


4 . Connections with the problem of directed polymers in 
random medium 


In the single particle case, the energy denominators 
associated to different sites along the paths are indepen¬ 
dent variables. Thus, the expression for the wave func¬ 
tion amplitude in FA, Eq. resembles the expression 
for the partition function of a directed polymer (DP) in 
a random potentiap2Ha4|^ 

with the thermal weights for 
the polymer configurations given by the amplitudes of 
the different paths. This analogy is not straightforward, 
due to the occurrence of negative contributions in (12|. 


Nevertheless, it has been fruitfully exploited both for the 
single particle problenPsHm^ and for problems of inter¬ 
acting spins on the Bethe lattic^SlHlIl 
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FIG. 6. (Color online) Probability of resonances 
P {Zr > \og{W/t)) for the variable Zr defined in Eq. (391 and 
d = 3. Asymptotically in r, the probability reaches 0 expo¬ 
nentially fast in the localized phase, and it reaches 1 double 
exponentially fast in the delocalized phase, in agreement with 
the analytic computations in Appendix In the plot, the 
squares are the results of the transfer mat rix calculation, the 
points of the dominating path (see Sec. |IV[ ) while the continu¬ 
ous lines are the exponential or double exponential fits. Very 
similar results are obtained for higher dimensionality. 


FIG. 7. (Color online) Power law divergence of the length 
scale 1{W) dehned in Eqs. ( [43| ) and ( |44| |. The values of 1{W) 
for fixed W are determined from hts such as the ones in Fig.[^ 
The power law fit produces a critical exponent iz ~ 1 and a 
critical value H4 compatible with the ones listed in Table [III 
The results shown here are for d = 3; very similar results are 
obtained for higher dimensionality. Notice how in the delo¬ 
calized phase the distance to observe a resonance is typically 
larger (for the same \W — Wc\) than the localization length 
in the localized phase. 


Motivated by this analogy ,the authors of Ref. 15^ have 
proposed a scaling form analogous to (40) for the log¬ 
arithm logg of the conductance of an Anderson model. 
There, the conductance in d = 2 is obtained from the 
Green functions, which are computed numerically within 
a modified FA, the modification consisting in taking en¬ 
ergy denominators that are not arbitrarily small but are 
bounded from below.l^ It is shown that the fluctuations 
of logg scale with an exponent uj{d = 2) = 1/3, and that 
the distribution of the variable u is compatible with a 
Tracy-Widom distribution. 

These results are consistent with the conjectur^^Ulll 
that in the strongly localized phase, where the expansion 
in non-repeating paths is best controlled, the Anderson 
model in dimension d belongs to the same universality 
class of the directed polymer in dimension 1 + D, with 
D = d — 1. In particular, the conjecture implies that 
in the limit of large r the distribution of logg has the 
scaling form (40), with uj{d) coinciding with the droplet 
exponent!^ in 1 -|- (d — 1) dimensions (which is exactly 
knowrP^to be equal to 1/3 for H = 1), and u having the 
same distribution of the fluctuations of the free energy in 
the disordered phase of the polymer (distr ibuted accord¬ 
ing to the Tracy-Widom distributiorF^^^ in D = 1). 


The values of the scaling exponents extracted from our 
data do not compare well with the droplet exponents 
uj{D = d — 1) of the DP, see Table |T] Curiously, they 
compare within errors with lo{D -1-1). We do not have 
an explanation for this curious behavior, and we leave its 
analysis for future work. Broadly speaking, the discrep¬ 
ancies with respect to the directed polymer results are 


generated by the fat-tail of the distribution of the paths 
amplitudes in (32), produced by the arbitrarily small en¬ 
ergy denominators. It might be that the finite size effect 
are more pronounced in the case of unbounded denomi¬ 
nators. On the other hand, it is quite natural to expect 
that the models of non-repeating paths with bounded 
amplitude considered in Ref. [3S] exhibit a stronger de¬ 
pendence on the dimensionality, due to the fact that the 
domination by one single path is less pronounced in that 
case. We comment more on this point in Sec. |IV[ 


C. Heisenberg model with random fields 

In order to test the forward approximation on a many- 
body problem, we consider an XXZ spin-1/2 chain in 
random magnetic field, 

L L L 

i—1 i—1 i—1 

(45) 

where periodic boundary conditions are assumed (s“ = 
s2-i-i)> ^'^'4 f4ie random fields hi are uniformly distributed 
in [—h,h]. This spin Hamil tonian (|4^ has been studied 
in a large number of work^^^SESHSST^ which numerical 
evidence of the existence of a localization-delocalization 
transition is provided, mainly based on exact diagonal- 
ization results. The critical disorder is estimatecP^ to be 
he — 3.72(6) for states in the middle of the energy band 
and parameters t = 1 and A = 1. 

As mentioned in the Introduction, the many body 
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FIG. 8. (Color online) Graph corresponding to the config¬ 
uration space of the XXZ spin chain, see Eq. (451, of total 
length L = Q and with periodic boundary conditions. Each 
site in the graph is associated to a product state in the basis 
of Si operators. Only sites corresponding to states with zero 
total spin are represented. The initial Neel state | 4,t • ■ •) ^-nd 
the final, totally flipped state | • ■ •) a-re highlighted with 

circles. The red, dashed edges form one of the shortest paths 
connecting the two states, of length L/2 = 3. 


problem can be seen as a single particle hopping problem 
in the “configuration space”. The latter is composed of 
the 2^ product states in the basis of sf, which span the 
full Hilbert space and diagonalize H (0). We denote these 
basis states with |n), and refer to them as the “compu¬ 
tational basis”. The mapping to an hopping problem is 
obtained by interpreting each state \n) as a vertex n of 
a graph, with associated random energy En defined by 
H{0)\n) = En\n). The third term in (45) provides the 
hopping between different sites, thus defining the geome¬ 
try of the graph. Note that due to spin conservation, the 
full configuration space, and consequently the graph, are 
partitioned into disjoint sectors corresponding to differ¬ 
ent total spin; we restrict to the sector of total spin equal 
to zero, corresponding to a connected graph with { 1 ^/ 2 ) 
vertices. See Fig. for a pictorial representation of the 
graph for L = 6. 

The effective hopping problem can be analyzed using 
the procedure set up in Sec. |II A| the amplitude dta of an 
eigenstate of the effective single particle problem is given 
in forward approximation by: 


In the many body language, Eq.(46) provides the ex¬ 
pression of the coefficients of the eigenstates of (45) in 
the computational basis, to lowest order in the coupling 
t. The exponential decay of the coefficients implies local¬ 
ization in the configuration space, meaning that the full 
many-body eigenstates are effectively a superposition of 
product states which differ only by configurations dis¬ 
tant 0(1) from the initial configuration. The two main 
consequences of this structure of the wave function is 
that they have si gnificantly less entanglement than er- 
godic stated^lMMI and that, using Kubo’s formula for lin¬ 
ear respons^i^, one can prove that they cannot support 
transport on macroscopic distances. 

Similarly to the Anderson case, we fix an initial config¬ 
uration of spins and we look at the amplitude in pertur¬ 
bation theory on the most distant, fully flipped configu¬ 
ration. In particular, we fix the localization center to be 
the site correspondent to the Neel state \ni) = | 4,t ...), 
and consider the wave function amplitude on the site cor¬ 
responding to the fully flipped Neel state |n 2 ) = | fi • ■ • 
These two sites ni, 712 are connected by 2(L/2)! paths on 
the graph, of length r = L/2 each. 

By means of the transfer matrix we compute the 
rescaled amplitude 


Zrih) = 


logl^.P 

2r 


(47) 


for different disorder strength h, with Ikj. given by ( [4^ . 
We consider spin chains of size 6 — 20 with hopping and 
interaction constants respectively t = 1 and A = 1, and 
= 1 — 6. Note that, despite the general framework is 
the same as in the Anderson problem, the transfer matrix 
calculation is by no means identical; indeed, in the many 
body case the energies associated to the different graph 
vertices are a linear combination of the independent ran¬ 
dom fields, and are thus correlated. Moreover, the num¬ 
ber of paths connecting two sites proliferates with the 
size of the chain L, with a scaling that is faster than 
exponential. These paths present correlations that are 
much stronger with respect to the Anderson problem, as 
we shall discuss in more detail in Sec. IIVI 


1. Distribution of the wave function amplitudes and critical 
disorder 

In Fig. I^we show the probability density of Zr{h) for 
a chain of length L = 20 and different values of h. Since 
in the many body case it is not possible to simplify the 
dependence on the disorder strength h, the criterion for 
the transition reads 

(Zoo(Lc)) =-logt, (48) 


*“("">= S n5;y(rE;. W 

pGspaths(ni,n2) nGp ^ 

where it is assumed that the eigenstate satisfies 4'„(n) —)■ 
Sn,ni for t 0. 


where {Zac,{h)) is the extrapolated value of the average 


of (47) for fixed h. Plots of {Zao{h)) are given in Fig. 10 


Here(Zoo(/i)) is extrapolated from the finite size values 
using the fitting function 


r {Zr{h)) = ci +{Zoo{h))r + C 2 r ^ 


(49) 
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FIG. 9. (Color online) Probability density of the random 
variable Zr{h) defined in Eq. (471, for an XXZ spin chain of 
length L — 20 (corresponding to r = 10) and different values 
of disorder h. Each curve is obtained with 3 • 10^ realizations. 



h 


FIG. 10. (Color online) Extrapolated value of the mean (Zoo) 
of the variable defined in Eq. (471. The crossing with 0 signals 
the many body localization/delocalization transition for t = 1 
(see Eq. (481). The error bars are obtained from the fitting 
procedure (see Inset). The resulting transition value is /ic = 
4.0T0.3. Inset. Finite size scaling of (r Zr) with the distance 
r between the Neel states ni and n 2 . The plot corresponds to 
h = 1. The fit is linear with an r~^ correction, see Eq. ( |49[ ), 
with parameters ci = —7.2 ± 0.4, (Zoo(2)) = 1-23 ± 0.02 and 
C 2 = 8.8 ± 0.7. The finite-r values for the mean are obtained 
over at least 10^ realizations for r < 7 and at least 2 • 10^ 
realizations for r > 7. 


the algorithm for the best path is not applicable in this 
context, and the limited system sizes accessible with the 
transfer matrix do not allow to investigate whether a 
scaling form exists also for (47) in the limit of large r. 
For the available system sizes, the distributions of the 
rescaled variables Zr{h) = {Zr{h) — {Zr{h))) jaz^th) do 
not seem to collapse to a unique curve, and the scaling 
of the variances (t| with r appears to be compatible 
with a power-law, but with exponent depending on the 
disorder strength h. However, a more refined numerical 
analysis is necessary to draw a conclusion on the asymp¬ 
totic behavior. 


2. Divergent length scales and critical exponents 


Fig.[n] shows the behavior of the probability of res¬ 
onances P {Zr{h) > — logt) as function of the distance 
between the Neel states. As expected, the r-dependence 
changes with the disorder: the probability decays to zero 
at large h, and increases towards one for the smaller h. 
We expect the convergence to be exponential in r on both 
sides of the transition; however, the exponential behavior 
is not clearly detectable in the delocalized phase, due to 
the few accessible system sizes. For h < we extract a 
length scale 1(h) by fitting the curves in Fig. 11 with the 
function: 

P (Zr(h) > 0) = a2(h) + H—(50) 

1 (h) r 

In the localized phase we perform instead the exponen¬ 
tial fit: 


P (Zr(h) > — logt) =ai(/i)exp 


1 (h)) ■ 


(51) 


The length scales 1(h) extracted with this procedure are 
shown in Fig. 12 together with the power law fit 1(h) = 
c\h — hc\~'^■ The fit is performed separately for h < he 
and h > he, resulting in an exponent close to 1 in both 
cases (see Fig. 12 for details). Note the asymmetry of the 
curve with respect to he, which indicates that at fixed 
\h — hc\ the typical distance to find a resonance in the 
delocalized phase is larger than the localization length at 
the corresponding value of disorder in the localized phase. 
A possible consequence of this phenomenon, which occurs 
also in the Anderson model (see Fig. [^, could be a large 
“critical region” in the dynamics in the delocalized phase. 


For t = 1, the critical point he is estimated from the 
condition (Zoo (he)) = 0. The resulting value is he = 
4.0 ± 0.3, which is, as expected, larger than the result 
derived with exact diagonalization. Notice also that the 
corrections oc are consistent with the intuition that 
the ojpA = 0 is the correct mean-field scaling for MBL. 


As we discuss in Sec. EYi in the many body case the 
sum (46) is no longer dominated by a single path. Thus, 


IV. THE STRUCTURE OF THE DOMINATING 
PATHS 


As mentioned in section Sec. ImBj the wave function 
amplitudes in FA can be interpreted as the partition func¬ 
tion for a directed polymer in random medium. However, 
a relevant difference is that while the weight associated 
to the polymer is bounded from abov^^, the single fac¬ 
tors in (33) are unbounded, with diverging average. As 
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In this section, we compare the statistics of the wave 
function amplitudes in FA with that of the optimal path, 
i.e. the path with maximal amplitude, both for the single 
particle problem and for the XXZ chain. We show that 
while in the first case the full sum is strongly dominated 
by the extremal path amplitude, in the many body case 
most of the paths have comparable amplitude. However, 
the much stronger correlation between them gives rise 
to non-negligible interference effects, resulting in strong 
cancellations. 


A. The single particle case 


FIG. 11. (Color online) Probability of resonances P{Zr(h) > 
— log t) as a function of the distance r between the two Neel 
states ni and n 2 , for t = 1. Asymptotically, the probability 
reaches zero exponentially in the localized phase and one in 
the delocalized phase. We average over 10^, 5 • 10® and 3 • 10® 
realizations for r < 8, r = 9 and r = 10, respectively; the 
plotted values of the disorder h are: h = 1 (points), h — 2 
(squares), h = 3 (diamonds), /i = 4 (upward triangle), h = 5 
(downward triangle) and h = 6 (circle). Linear and exponen¬ 
tial fits in the delocalized and localized regions respectively 
are plotted as continuous lines, see Eqs. (501 and (511. 


500 



2 3 4 5 6 


h 

FIG. 12. Color online) Divergence of the length scales l{h) 
extracted from the fits of the probability of resonances as a 
function of r. The vertical dashed line indicates the critical 
value he obtained in Fig. The dotted curve is a power 
law fit of the form c\h — hc.\~'', resulting in a critical exponent 
i/L ~ 1.12 ± 0.06 for h < he and vr = 1.1 ± 0.2 for h > he- 


a consequence, the strongly localized phase in the An¬ 
derson model (where the FA is best controlled) always 
corresponds to a “frozen” phase of the directed polymer, 
in which most of the weight in the total sum (12) is given 
by one single path. An interesting questiorP^ is whether 
the freezing phenomenon persists in the delocalized phase 
in some form, and what this implies t hat th e structure 
of the eigenstates close to the transitiorP^^^. 


For the finite dimensional case, we compute the ampli¬ 
tude oj* of the optimal path p* dominating the sum (12) 
by means of the Dijkstra algorithnP^, a graph-search al¬ 
gorithm that determines the path minimizing a given cost 
function. We consider a directed cube (such as the one 
in Figure and assign a positive cost y to each directed 
edge (ij): 


X{i,j) = log|ej| -min{log|efc|}. 

k 


(52) 


The total cost of a path p is the sum of the costs of 
the edges belonging to it, and the path p* with maximal 
amplitude is the one minimizing the total cost function. 
In order to compare with the transfer matrix results, we 
compute the ratio between uj* and the full sum (32) com¬ 


puted via the transfer matrix technique, for the same 
given disorder realization. The distribution of the ratios 
turns out to be very narrowly peaked around one. Fig¬ 
ure displays its average as a function of the length of 
the paths r for d = 3, which is extremely close to one, 
uniformly in the path length. 

As a further check of the agreement between the val¬ 
ues computed with the two methods, we plot in Fig. 
the r-dependence of the probability © with d = 0, de¬ 
termined with the substitution |'0j.| —> uj*. The data 
are plotted as points, which are almost indistinguishable 
from the transfer matrix results (squares). This indicates 
that the statistics of distant resonances is fully captured 
by the optimal path. Thus, in the single particle case the 
correlation between different paths does not play a rele¬ 
vant role, since the sum is dominated by the extremum, 
as it would happen for independent random variables 
with fat-tailed distribution. Based on this observation, 
the numerical analysis outlined in the previous sections 
can be carried out for much bigger system sizes with re¬ 
spect to the ones accessible with the transfer matrix tech¬ 
nique, since the Dijkstra algorithm has lower complexity 
than the transfer matrix (indeed, the time complexity is 
~ 0{e + ulogu) and the space complexity is ~ O(u^), 
where v is the number of vertices and e is the number of 
graph edges). 

We perform the same analysis also for the modified 
forward approximation discussed in Ref. 1361 by tak¬ 
ing the energy denominators uniformly distributed in 
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FIG. 13. Average ratio between the dominating path weight 
uj* computed with the Dijkstra algorithm, and the sum in Eq. 
( |32[ | computed using a transfer matrix technique. The ratio 
is taken between values corresponding to the same disorder 
realization. The plot corresponds to d = 3, and each point is 
averaged over 3 • 10"^ disorder realizations. Similar results are 
obtained for higher dimensionality, for those r accessible with 
the transfer matrix technique. The standard deviation error 
bars are within the point size. 


[—1,—FF“^] U [TT“^,1] in d = 3 for two values of the 
cutoff, FT = 25 and FF = 35. We find that in this case the 
ratio between the maximal path and the transfer matrix 
result departs from one for increasing r. This suggests 
that more than one path dominates the transfer matrix 
result. It is natural to expect that in this case the num¬ 
ber of dominating paths depends on the geometry of the 
system, thus introducing a stronger dependence on the 
dimensionality, see also the comments in Sec. IIIB| 

For the case of unbounded denominators, we compute 
the inverse participation ratio (IPR) of the edge weights 
contributing to w*, for G [—1,1] (i-e. FF = 2). We 
define 


IPR = 


(E.ioglGl)^ 

E.(loglGl)^ 


(53) 


where i labels the sites belonging to the optimal path p*. 
We find that the disorder-averaged IPR scales linearly 
with the length of the path r, indicating that an exten¬ 
sive (in r) number of edges contributes to the total path 
weight, and cooperate to produce the atypically big path 
weights dominating (32). Fig. 14 shows the distribution 


of the absolute value of the energies along the optimal 
path for FF = 2, d = 3, r = 210 and = 0. The fitting 
function has the form 


Pr{e) = Cr + . (54) 

The power-law behavior is consistent with the consider¬ 
ations in Ref. |42l Adapting their reasoning to the finite 
dimensional case, one can argue that asymptotically in r 
(and under the hypothesis of independent paths) the bi¬ 
ased energy distribution along the optimal path has the 


FIG. 14. (Color online). Probability distribution p(e) of the 
energy denominators along the optimal path, see Eq. (55l. 
The plot corresponds to d = 3 and r — 210. The dashed red 
line is the fitting function of Eq. (541, with fitting parameters 
= -0.95 ± 0.04, br = 1.04 ± 0.03 and = 0.472 ± 0.005. 
Very similar results are obtained for higher dimensionality. 
Inset. Plot of the exponents in the fitting function of 
Eq. (54), as a function of r. Due to the absence of a theoret¬ 
ical reasoning for the finite size scaling, we fit the curve con¬ 
sidering logarithmic and 1 / y/r corrections. The green small- 
dashed curve is a fitting function of the form a+c/ log(r), with 
fit parameters a = —0.73 ± 0.05 and c = —1.4 ± 0.3; the red 
large-dashed curve is a fitting function of the form a + cj y/r, 
with fit parameters a = —0.57 ± 0.02 and c = —1.4 ± 0.3. 
The asymptotic value a obtained with the logarithmic fitting 
function is compatible with the solution of the equation (56) 
for d = 3. 


form 


P(e) = 


1 - 2a; 


(55) 


with X solving the d-dependent equation 


log 


l-2a: 


2a; 


1 - 2a; 


= 0 . 


(56) 


Fitting the r-dependence of the coefficients Cr,6r,«r one 
finds that the asymptotic limits are in agreement with 


(55); for details see the inset of Fig. 14 


B. The many-body case 


When performing the same type of analysis for the 
Heisenberg chain, we find that the statistics of the sum 
(46) is not well reproduced by the optimal path alone: 
the distribution of the ratios between the full sum and 
the optimal path is very wide and peaked at values that 
are far from one. Thus, despite also in this case the 
amplitude of the single paths are fat-tailed distributed, 
there is not a single one dominating. Instead, we find 
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FIG. 15. Average IPR* of the paths, Eq. ( |57[ ), as a function 
of the number of paths N* = 2(L/2)! for random spin chains 
of different lengths L. The IPR* is linear in the total number 
of paths. 


that the average IPR* of the paths amplitudes (which 
we denote with Wr,): 


IPR* = 


(s. 


E 


P 


(57) 


scales linearly with the total number of paths N* = 
2{L/2)\, indicating that there are factorially many (in 
the length of the chain L) paths having amplitudes that 
are comparable in absolute value. This is a signature of 
the strong correlations between the paths, which is not 
surprising in view of the many-body nature of the model. 
Following Ref. [HI one can argue that the strongest cor¬ 
relations are among those paths associated to processes 
in which the same spin flips occur, but in different or¬ 
der: the different orderings of the flips produce different 
energy denominators in (46), and thus different path am¬ 


plitudes; however, the resulting terms are correlated, and 
one can expect that for those realizations of the random 
fields producing one particularly large path weight, the 
other ones (related to it by permutation of the order of 
the spin flips) will also have a large amplitude in abso¬ 
lute value. However, in the sum (46) the paths contribute 


with well defined relative signs, leading to cancellations 
between these factorially many terms (see Ref. [T6| for an 
explicit calculation for a model of interacting fermions), 
which are fully taken into account only with the transfer 
matrix method. 


V. CONCLUSION 

In this work we have discussed the advantages and the 
limitations of the forward approximation applied to both 
single- and many-body quantum disordered systems. In 


particular, the FA has been used to obtain an expression 
for the wave functions, that can be computed by means 
of a transfer matrix technique. The statistical analysis of 
the wave functions allows to determine the critical values 
of the disorder (exact in large d), the critical exponents 
of the localization length (which turn out to be mean- 
field) and the universal distribution of the eigenfunctions’ 
coefficients. 

For the single particle case, the amplitudes of the wave 
functions in FA turned out to be very well approximated 
by only one path, the dominating path: this has been 
exploited to investigate larger system sizes with respect 
to the ones accessible with the transfer matrix (the algo¬ 
rithm that computes the best path runs in time linear in 
the number of edges of the underlying graph, and is not 
as memory-demanding as either the transfer matrix of 
shift-invert exact diagonalization). The extremely good 
agreement (within statistical error) of the predicted crit¬ 
ical values for the disorder for d > 5 suggests that the 
approximation should be predictive also for the proper¬ 
ties of the wave functions at the critical point. We have 
not investigated in details these implications, but we fore¬ 
see the wave functions to have a sparse structure which 
is similar to that discussed in high-coordination Bethe 
lattice^^. The strong similarities with the problem of 
directed polymers in random medium have also been ad¬ 
dressed; however, from the statistical analysis it emerges 
that the scaling exponents describing the fluctuations of 
the wave functions are non-mean field, but also not equal 
to those of the directed polymer. Moreover, the limiting 
distribution of the appropriately rescaled wave functions 
seems to depend weaker on the dimensionality with re¬ 
spect to the directed polymer case. 

In many-body problem there is no concentration of am¬ 
plitude on a small number of paths, but there are strong 
cancellations between them: as a result, the full sum 
over factorially-many paths is only exponentially large 
(or small) in the system size. The correlation between 
the paths has been discussed in detail in Ref. |H1 and 
this work can be interpreted as a numerical test of the 
claims in that work. For the XXZ chain with random 
fields, the critical value predicted within the approxima¬ 
tion {he = 4.0 ±0.3) is very close to the most updated re¬ 
sult obtained with exact diagonalization; thus, this tool 
furnishes an alternative route to exact diagonalization, 
that can be applied to significantly larger system sizes. 
We leave to future work the question of how to incorpo¬ 
rate higher-order corrections in the FA within the trans¬ 
fer matrix scheme, and how they affect the critical expo¬ 
nents, and the accuracy of the results. 

As a conclusive remark, we would like to briefly com¬ 
ment on the nature of the FA as a mean-field approxima¬ 
tion for Anderson localization. For certain, the value of 
the critical disorder Wc for the onset of localization grows 
indefinitely with d, and in high d the hopping t becomes 
an almost negligible perturbation at the transition. The 
fact that the error in Wc essentially disappears around 
d ~ 6 is a strong indication of this. This feature is quite 
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peculiar, since in ordinary, second order phase transitions 
the critical exponents above the upper critical dimension 
are correctly reproduced by the mean held approxima¬ 
tion, but the location of the transition (e.g. the critical 
temperature) is not. In this sense the locator expansion 
(that to lowest order reduces to the FA considered in this 
work) becomes a better suited candidate for a mean held 
than the 2 + e expansion of the nonlinear supersymmet¬ 
ric sigma model (NLScrM))^. It is plausible that there 
is a held-theoretical description of the FA which can be 
put in direct relation with the NLSctM; this is an obvi¬ 
ous direction in which to continue this work. In addition, 
the relation with the Bethe lattice results can be further 
investigated, given that the critical Wc predicted by the 
FA does not correspond to that of a Bethe lattice of any 
(integer) coordination number. 
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Appendix A: Probability density of Zr 

An estimate for the probability density of Zr can be 
obtained making use of the fact that the sum over (33) 
is dominated by the single path with maximal weight: 


max 

pGspaths 


log 


(Al) 


iGp 


If the correlations between the different path weights 
are neglected, the calculation is similar to the one per¬ 
formed for the Bethe lattice case. In particular, one hnds 
for the cumulative function of Zr the following expres¬ 
sion: 


P {Zr < a) = exp 


V^logfl- ^ [ e *dt] 

\ {r -V'- Jr{a-log2) J 

(A2)- 

where Nr ~ (F is the total number of paths on which the 
maximum is taken. This implies the following form for 
the probability density of Z'r = Zr — log 2: 

(A3) 

where we introduced the monotone decreasing function 


Iriz') = 


(r- 1)! 


P *dt. 


(A4) 


The typical value of Z', denoted 2 ;*, is dehned by the 
equation 


Nrlriz*) = Nr 


y^T— 1 


(r- 1)! 


t 


r— 1 ^—rt 


dt = 1. 


(A5) 


The solutions of ( |A5| ) approach a finite limit z* for 
r —> 00 , which is related to the critical value of disorder 
by 2 * = log(IFc/(2t)), as previously discus sed. Using 
that Nr ~ dr and computing the integral in (A5) with a 
saddle point calculation (assuming z* > 1), 


one recovers 


the condition (28) for Wc, with the substitution K ^ d. 


For increasing r the probability density of Z{. peaks 
at the typical value, with tails that approach zero in the 
limit r —> 00 . In particular, for z' > z* the decay of 
the tail is exponential in r. Indeed, in this regime the 
product Nr Ir{z) is itself exponenti ally decreasing with r; 
thus, the rightmost exponential in (A3) rapidly converges 


to one, and the distribution Pr{z') approaches zero with 
a tail of the form 


Pr{z') ~ Q-'^{^'-^°Sz'-loeide))+o{r) 


(A6) 


When z' becomes smaller than the typical value z*, the 
product Nr Ir{z') increases exponentially. Since for large 
r the integral Ir{z') is still exponentially small for all 
z' > 1 -I- 0(l/r), one can still set log [1 —/^(z:')] ^ 
—Ir{z') ~ exp (—rz' -I-r log(ez')-|-o(r)). Thus, in this 
regime the probability density of Z{. decays to zero much 
faster, double-exponentially with r 

Pr{z') ~ exp (-dJ' -t 0(r)) . (A7) 


Note that (for r large enough) the interval in which 
1 < z' < z* does not shrink to zero for d > 3 , giv en 
that the value z* obtained from the condition (A5) is 
always bigger than one. When z' approaches one, the 


probability in (A4) is no longer a large deviation prob¬ 


ability, i.e. it is no longer exponentially small in r: the 
term log [1 — Ir{z')\ approaches a constant function of z', 
and the main scaling is given by the factor dT. 

Finally, exactly at z' = z*, using (A5) and performing 


the integral with an integration by parts, one finds that 
the probability density can be written as 


Pr{Zr) = 


1 - d-'' 
exp I —1 — 


1 + E 


(r- 1)! 


71 = 1 
1 1 
^ ~ 3dW 


(r — 1 — n)!r’ 


:«)■ 


-1 


(AS) 


which diverges like r when r —>■ 00 . 

Given the tails of the distribution of Zr computed in 
this approximation, it is immediate to derive the asymp¬ 
totic decay of the probability of resonances in the local¬ 
ized phase. Indeed, for W > Wc, the prob abili ty (35) is 
is a large deviation for Zr ■ Making use of (A6) we find 


p(Zr> log = r ^-r{^'-^og.'-\og{de))+o{r)^^, 

\ t J Aog(D) 

(A9) 
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with 


1 , ( W 1 \ 

1{W) ^^\2tdelog{W/2t)J ' 


(AlO) 


Thus, within this approximation for W approaching 
Wc from above one finds 


1{W) 


W-W^’ 


(All) 


thus the length scale diverge s at the transition with a 


critical exponent equal to 1. 

Simil arly, for 2te < IT < Wc, making use of (A2) and 
of (A7| we find: 


P 


(^Zr < log 



« exp 


2 tde 

~W~ 


log 




= exp + 0(r)) . 

(A12) 
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